R/sink fit.R

In mlrioux/sinkworx: Import, Model, and Plot Cell Sinking Data

# June 4, 2019

#'@title Run segmented regressions on cell sinking data
#'@description Fits a segmented regression to cell sinking data (square root RFU values) by well. The model parameters can then be used to calculate sinking rate.
#'@usage sink_fit(database_rfu, breakpoints = c(20, 50), attempts = 19, export_output = FALSE,
#'save_as = "Sinking_Fits.csv")
#'@param database_rfu dataframe containing the columns "Plate", "Well", "Elapsed.Time.m", and "RFU". Can be generated by the user, or by running the import functions on raw data files.
#'@param breakpoints initial x-axis estimates where the regression should be broken up.
#'@param attempts number of times the segmented regression should be run until a fit is generated. This is for internal error handling. Do not change unless needed.
#'@param export_output whether the resulting dataframe should be exported to .csv.
#'@param save_as desired filename, if results are to be exported.
#'@return A dataframe containing the model parameters, with each row representing a single well. The output parameters include well bottom (mean of 3 lowest points after second breakpoint), slopes 1, 2, and 3, breakpoints 1 and 2, y-intercepts on each of the slopes, and 95% upper and lower confidence intervals on each slope.
#'@examples
#'sink_fit(database_rfu, breakpoints = c(2, 10), export_output = T, save_as = "Sink_Parameters.csv")
#'@note The parameters generated by this function are NOT final sinking rates. You may use the sink_rate() function to generate sinking rate estimates from the parameters, or you may wish to apply your own equation.

# ------------------------------------------------------------------------------

sink_fit <- function(database_rfu,  breakpoints = c(20, 50), attempts = 19,
                     export_output = FALSE, save_as = "Sinking_Fits.csv") {

# ------------------------------------------------------------------------------

  library(magrittr)

  plot_by <- paste(database_rfu$Plate, database_rfu$Well)

  wells <- unique(plot_by)

  Sink_Fits <- NULL

  for (w in 1:length(wells)) {

    focus_well <- wells[w]

    well_data <- dplyr::filter(database_rfu, plot_by == focus_well)

    # Start with linear model

    sink_lm <- lm(sqrt(RFU) ~ Elapsed.Time.m, data = well_data)

    # Now attempt segmented regression until it is no longer null

    sink_seg <- NULL

    attempt <- 0

    while (is.null(sink_seg) && attempt <= attempts) {

      attempt <- attempt + 1

      try(sink_seg <- segmented::segmented(sink_lm, seg.Z = ~ Elapsed.Time.m, psi = list(Elapsed.Time.m = c(breakpoints))))

    } # close attempt loop

    slope_1 <- coef(sink_seg)[[2]]

    slope_2 <- coef(sink_seg)[[2]] + coef(sink_seg)[[3]]

    slopes <- broom::tidy(segmented::slope(sink_seg)$Elapsed.Time.m)

    slope_3 <- slopes$Est.[[3]]

    breakpoint_1 <- sink_seg$psi[[3]]

    breakpoint_2 <- sink_seg$psi[[4]]

    intercept_1 <- coef(sink_seg)[[1]]

    intercept_2 <- (intercept_1 + slope_1 * breakpoint_1) - (slope_2 * breakpoint_1)

    intercept_3 <- (intercept_2 + slope_2 * breakpoint_2) - (slope_3 * breakpoint_2)

    slope_1_LCI <- slopes$CI.95...l[[1]]

    slope_1_UCI <- slopes$CI.95...u[[1]]

    slope_2_LCI <- slopes$CI.95...l[[2]]

    slope_2_UCI <- slopes$CI.95...u[[2]]

    slope_3_LCI <- slopes$CI.95...l[[3]]

    slope_3_UCI <- slopes$CI.95...u[[3]]

    bottom_RFU <- dplyr::filter(well_data, Elapsed.Time.m > breakpoint_2) %>% dplyr::arrange(sqrt(RFU))

    if (nrow(bottom_RFU) >= 3) {

      well_bottom <- mean(sqrt(bottom_RFU$RFU)[0:3])

    } else if (nrow(bottom_RFU) < 3) {

      well_bottom <- mean(sqrt(bottom_RFU$RFU)[0:2])

    } # close well bottom calculation

    sink_fits <- as.data.frame(cbind(well_data[1, ],
                                     well_bottom,
                                     slope_1, slope_2, slope_3,
                                     breakpoint_1, breakpoint_2,
                                     intercept_1, intercept_2, intercept_3,
                                     slope_1_LCI, slope_1_UCI,
                                     slope_2_LCI, slope_2_UCI,
                                     slope_3_LCI, slope_3_UCI)) %>%
      dplyr::mutate(., Elapsed.Time.m = NULL, RFU = NULL)

    Sink_Fits <- rbind(Sink_Fits, sink_fits)

  } # close loop through wells

  if (export_output == TRUE) {

    write.csv(Sink_Fits, file = save_as, row.names = F)

  } # close export option

  return(Sink_Fits)

} # end of function